﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.EisPack
{
    /// <summary>
    /// Forms the eigenvectors of a real general matrix by back transforming those of the corresponding 
    /// upper hessenberg matrix determined by elmhes.
    /// </summary>
    [Serializable]
    public static class ElmbakClass
    {
        /// <summary>
        /// Forms the eigenvectors of a real general matrix by back transforming those of the corresponding 
        /// upper hessenberg matrix determined by elmhes.
        /// </summary>
        /// <param name="lo">The lo parameter.</param>
        /// <param name="hi">The hi parameter.</param>
        /// <param name="a">contains the multipliers which were used in the reduction by elmhes.</param>
        /// <param name="intch">Contains information on the rows and columns interchanged in the reduction by elmhes.</param>
        /// <param name="m">The number of columns of z to be back transformed.</param>
        /// <param name="z">Contains the real and imaginary parts of the eigenvectors to be back transformed in 
        /// its first m columns.</param>
        public static void Elmbak(int lo, int hi, double[,] a, int[] intch, int m, double[,] z)
        {
            /* This routine is a translation of the Algol procedure from
             * Handbook for Automatic Computation, vol. II, Linear Algebra,
             * by Wilkinson and Reinsch, Springer-Verlag.
             */
            double x;
            int i, j, kp1, la, mm, mp, mp1;

            if (m == 0)
            {
                goto _200;
            }
            la = hi - 1;
            kp1 = lo + 1;
            if (la < kp1)
            {
                goto _200;
            }
            for (mm = kp1; mm <= hi; mm++)
            {
                mp = lo + hi - mm;
                mp1 = mp + 1;
                for (i = mp1; i <= hi; i++)
                {
                    x = a[i, mp - 1];
                    if (x == 0.0)
                    {
                        goto _110;
                    }
                    for (j = 0; j <= m; j++)
                    {
                        z[i, j] += (x * z[mp, j]);
                    }
                    _110:
                    ;
                }
                i = intch[mp];
                if (i == mp)
                {
                    goto _140;
                }
                for (j = 0; j <= m; j++)
                {
                    x = z[i, j];
                    z[i, j] = z[mp, j];
                    z[mp, j] = x;
                }
                _140:
                ;
            }
            _200:
            ;
        }
    }
}